RD-R1<4  12S  REVIEH  IMPLEMENTATION  AND  TEST  OF  THE  QAZID 

COMPUTATIONAL  METHOD  NITN  R  VIEH  TO  HAVE  ROTOR 
APPLICATIONSCU)  NAVAL  POSTGRADUATE  SCHOOL  HONTEREV  CA 
UNCLASSIFIED  T  F  SALACKR  DEC  83  NPS<7-85-ei2  F7G  2074 


HTIC  FILE  COPY  AD-A164  125 


NPS67-85-0 

NAVAL 


POSTGRADUATE  SCHOOL 

Monterey,  California 


DTIC 

f-TLECTE 
FEB  1  4  1986  ^ 

B 


THESIS 


REVIEW, 

QAZID 

VIEW 

IMPLEMENTATION  AND  TEST  OF  THF 
COMPUTATIONAL  METHOD  WITH  A 

TO  WAVE  ROTOR  APPLICATIONS 

by 

Thomas  F. 

Salacka 

December 

1985  . 

Thesis  Advisor; 

Ray  P.  Shreeve 

Ap:.roved  for  public  release;  distribution  is  unlimited 


Prepared  for:  Naval  Air  Systems  Command 
Washington,  DC  20361 


211  d2Q 


NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


Rear  Admiral  Robert  H.  Shumaker  David  A.  Schrady 

Superintendent  Provost 


This  thesis  prepared  in  conjunction  with  research  sponsored 
in  part  by  the  Naval  Air  Systems  Command,  Washington,  DC,  under 
Air  Task  #A310310E/186A/WR024-03-001 . 

Reproduction  of  all  or  part  of  this  report  is  authorized. 


Released  By; 


UNCLASSIFIED 

SEC'jRirv  Ct.ASSi^'CATI(>N  6^  fwiS 


/!/)/) 


la  REPORT  SECURITY  CLASSIFICATION 

Unolassif i&d _ 

2a  SECURITY  CLASSIFICATION  AUTHORITY 
)  ;b  OEC'-ASSiFiCATiON/ DOWNGRADING  SCHEDULE 

I _ 

4  PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 

NPS67-85-012 


REPORT  DOCUMENTATION  PAGE 

jib.  RESTRICTIVE  MARKINGS 


3  DISTRIBUTION /Ay/AILA8ILITY  OF  REPORT 

Approved  for  public  release; 
distribution  is  unlimi ted 

S  MONITORING  ORGANIZATION  REPORT  NU'V13ER(S) 

NPS67-85-012 


I  6a.  .NAME  OP  PERFORMING  ORGANIZATION  6b  OFFICE  SYMBOL  7a.  NAME  OF  MONITORING  ORGANIZATION 
i  (If  ippikabi*) 

Naval  Postgraduate  School  Code  67  Naval  Postgraduate  Scl 


6c  ADDRESS  (Ofy,  state,  and  IlPCoda) 

Monterey,  California  93943-5100 


Naval  Postgraduate  School 

7b.  ADDRESS  (C/ty,  State,  and  ZIP  Code) 

Monterey,  California  93943-5100 


!  3a  NAME  OF  FUNDING /SPONSORING  Bb  OFFICE  SYMBOL  9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 

ORGANIZATION  (It  applicable) 

Naval  Air  Systems  Camand  Air  310E  N0001985WR  51147 

3c  ADDRESS  (C/ty,  State,  and  ZIP  Code)  10  SOURCE  OF  FUNDING  NUMBERS 

PROGRAM  I  PROJECT  I  yaS<  IwORK  JN 

Washington,  D.C.  20361  element  no  no  no  accessioi 

_ 61153N  WR024 _ 03 _ 001 

title  (Include  Security  Clatsification) 

REVIEW,  IMPLEMENTATION  AND  TEST  OF  THE  QAZID  COMPUTATIONAL  METHOD  WITH  A 
VIEW  TO  WAVE  ROTOR  APPLICATIONS _ 

■2  personal  AUTMOR(S) 

Salacka,  Thomas  F. 

■?J  TVPC  nc  report  i;b  TIME  COVERED  14  DATE  OF  REPORT  (Year, /Vfontft,  Day)  IS  P/iCiE  COlNT 

J-laster'^s  Thesis  from__ _ tq  1985,  December  119 

I  '6  Supplementary  notation  ^ 


PROJECT 

task 

NO 

NO 

WR024 

03 

work  unit 
ACCESSION  NO 


COSATI  COOES  I 

-  EcD 

GROUP 

SUB-GROUP 

18  SUBJECT  TERMS  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 


Dynamics  ;  Euler  /-  / 


9  abstract  {Continue  on  reverse  if  necessary  and  identify  by  block  number) 

The  QAZID  method  is  redeveloped  in  detail  and  implemented  in  a 

first  order,  one-dimensional  FORTRAN  program,  EULER-1.  The  program 

is  tested  on  the  shock  tube  problem  and  results  are  presented  for 

various  computational  meshes  and  initial  conditions.  Based  on  good 

results  of  the  EULER-1  code,  recommendations  are  made  for  future- — 


■  '“'3  J^IQN,  AVAILABILiTv  OF  abstract  21  ABSTRACT  SECTURITY  CLASSIFICATION 

g  .•.c.assif'eo.'jnl  mited  □  same  as  rpt  Ddtic  users  Unclassified 

i  .J  -.AME  OF  RESPONSIBLE  NDiViDUAL  22b  TELEPHONE  f/nc/ud*  Area  Code)  I  22c  OFFICE  SYMBOL 

JProf^^Ray^^^^^Shreeve^^^^^^^^^^^^^^^  (4£8)646|-££2^^i_____i__£2£^,_^2££«« 

00  FORM  1473,  84  MAR  33  APR  edition  may  be  used  until  exhausted  SECURITY  CLASSIFICATION  OF  this  PAGE 

All  other  editions  are  obsolete 

I  UNCLASSIFIED 


Approved  for  public  release;  distribution  is  unlimited 


Review,  Implementation  and  Test  of  the 
QAZID  Computational  Method  with  a 
View  to  Wave  Rotor  Applications 

by 

Thomas  Francis  Salacka 
Lieutenant,  United  States  Navy 
B.S.,  United  States  Naval  Academy,  1977 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


MASTER  OF  SCIENCE  IN  AERONAUTICAL  ENGINEERING 


from  the 

NAVAL  POSTGRADUATE  SCHOOL 
December  1985 


Author : 


Thomas  Francis  Salacka 


Approved  by: 


C  I 


Raymc 


ireeve. 


ivisor 


Max  Platzer,  Chairman, 
Department  of  Aeronautics 


Dean  of  Science  and  Enaineering 


ABSTRACT 


s' 


The  QAZID  method  is  redeveloped  in  detail  and  implemented 
in  a  first  order,  one-dimensional  FORTRAN  program,  EULER-1. 
The  program  is  tested  on  the  shock  tube  problem  and  results 
are  presented  for  various  computational  meshes  and  initial 
conditions.  Based  on  good  results  of  the  EULER- 1  code, 
recommendations  are  made  for  future  extensions  and  testing 
to  validate  the  suitability  of  the  QAZID  method  for  wave 
rotor  applications. 
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SYMBOLS  USED  IN  THE  TEXT 
Definition 
Speed  of  sound 
Static  pressure 

.Modified  Riemann  variable  (Q  =  q  +AS) 
Absolute  velocity  magnitude 
Reversible  heat  transferred 
Modified  Riemann  variable  (R  =  q  -AS) 

Gas  constant 

A  modified  form  of  entropy 

Unit  vectors  in  the  s,  n,  and  m  directions 

Static  Temperature 

Time 
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Vector  of  the  principle  variables,  Q,R,S 
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DARRAY 

DEL**H 

DELX (K) 

DET,**L 
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governing  equations 

LMD(K)  The  characteristic  trajectories  in  the 
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Small  spatial  change 
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Flow  angles  with  respect  to  reference 
coordinate  planes 
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The  average  value  of  A  between  the  limits  of 
integration  of  Z (K) 
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The  distance  from  the  point  where  the  K 
characteristic  crosses  the  Icnown  time  level  to 
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The  change  in  the  variable  **  in  going  from  node 
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J 

JSTOP 
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**  over  one  time  step 

The  change  in  the  velocity  magnitude  across  a 
shock,  (high-low) 

The  ratio  of  the  density  across  a  shock,  (low/high) 

A  value  which  indicates  the  precision  to  which  the 
characteristics  are  to  be  calculated 

The  computed  error  in  the  calculation  of  the 
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Parameter  which  controls  the  type  of  output  pro¬ 
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pressure  and  density,  2  =  Comparison  of  density 
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2/(G-l) 
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Subscript  which  denotes  the  spatial  node 

The  interpolated  value  of  **  at  the  point  on  the 
known  time  level  where  the  characteristic 

crosses 

The  result  of  integrating  Z(K) 

The  node  to  the  right  of  discontinuity  L 
The  time  level 

The  number  of  time  levels  desired  to  be  computed 

A  subscript  which  indicates  the  characteristic 
being  dealt  with.  1  =  Q+A  2  =  Q-A  3  =  Q 


K 


L  A  subscript  which  denotes  the  type  of  discontinuity. 

1  =  schock  2  =  contact  discontinuity  3  =  head 
of  rarefraction  wave  4  =  tail  or  rarefraction 
wave 

M  Parameter  which  controls  the  types  of  discontinui¬ 

ties  which  EULER-1  tracks.  1  =  shocks  only. 

2  =  shocks  and  contact  discontinuities 

The  number  of  nodes  in  the  computational  grid 

A  temporary  storage  location  for  variables  ** 
updated  in  time 

The  array  of  densities  plotted  by  the  graphics 
routines 

PR  The  ratio  of  pressures  across  a  shock,  (high/low) 

PRI  The  initial  pressure  ratio  across  the  diaphragm 

(high/low) 

*PRIM(K)  Suffix  which  indicates  the  spatial  derivative  of 
*  at  the  present  time  level,  where  *  is  either  A 
or  Q 

QARRAY  The  array  of  velocities  to  be  plotted  by  the 

graphics  routines 

QLI  The  initial  value  of  the  velocity  to  the  left  of 

the  diaphragm 

QRI  The  initial  value  of  the  velocity  to  the  right  of 

the  diaphragm 

QQJO  The  measured  change  in  the  Riemann  variable  QQ 

across  a  shock  (high-low) 

QQJE  The  change  in  the  Riemann  variable  QQ  across  a 

shock  calculated  analytically  (high-low) 

SIGMA (L,J)  Discontinuity  locations.  J  =  1  indicates  the 

known  time  level  and  J  =  2  indicates  the  unknown 
time  level 

The  change  in  the  primary  variables  with  respect 
to  time  during  one  time  step  (**  =  QQ,  RR,  S) 


N 

NEW** 

PARRAY 


**STEP 


SARRAY 

SKIP 

TRI 

VCDE 

VHEAD 

VTAIL 

VS 

VSE 

X 

XARRAY 

XINIT 


The  array  of  entropy  values  to  be  plotted  by  the 
graphics  routines 

The  number  of  time  steps  between  outputs 

The  initial  temperature  ratio  across  the  diaphragm 
(high/low) 

The  exact  velocity  of  the  contact  discontinuity 

The  exact  velocity  of  the  head  of  the  rarefraction 
wave 

The  exact  velocity  of  the  tail  of  the  rarefraction 
wave 

The  shock  wave  velocity  computed  by  EULER-1 

The  exact  velocity  of  the  shock  wave 

The  non-dimensional  spatial  position 

The  array  of  node  locations  used  by  the  graphics 
routines 

The  initial  location  of  the  discontinuities.  Used 
to  calculate  the  exact  solution 


X2  (L) 


The  position  of  the  node  to  the  right  of  a 
discontinuity 
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INTRODUCTION 


I  . 

A  program  is  underway  at  the  Naval  Postgraduate  School's 
(NPS)  Turbopropulsion  Laboratory  to  evaluate  the  wave  rotor 
concept.  The  wave  rotor,  operating  as  a  component  in  a  gas 
turbine  engine,  uses  unsteady  wave  propagation  in  tube-like 
passages  to  compress  incoming  air  before  it  goes  to  the  com¬ 
bustion  chamber.  The  combustion  chamber  output  is  routed  back 
to  the  rotor  to  create  the  unsteady  waves.  Thus  the  rather 
simple  rotor,  with  "partial  admission"  inlet  and  outlet  ports, 
acts  both  as  a  compressor  and  as  a  turbine.  The  alternation 
of  hot  and  cool  gases  through  the  same  hardware  aids  cooling 
and  allows  higher  operating  (cycle)  temperatures  to  be  used. 
The  interaction  of  the  hot  and.  cool  gases  within  the  passages 
of  the  wave  rotor  through  the  wave  patterns  that  are  created 
pose  a  difficult  problem.  The  work  underway  aims  to  develop 
and  validate  preliminary  design  and  performance  analysis  tools 

One  of  the  tools  needed  is  a  computer  program  which  can 
be  used  to  construct  wave  interactions  in  a  design  process 
and  then  accurately  predict  the  performance  of  the  design. 
Until  efficient  and  accurate  methods  of  solving  the  full 
Navier-Stokes  equations  for  unsteady  turbulent  flows  are 
developed,  the  design  program  will  be  based  initially  on  the 
solution  of  the  unsteady  Euler  equations  and  a  method  devised 
to  represent  losses.  Once  the  program  is  developed,  it  must 


be  verified  against  experiment.  This  is  the  overall  goal  of 
the  present  program  in  which  a  v’ave  rotor  apparatus  has  been 
assembled  and  methods  of  measuring  the  unsteady  pressures  and 
temperatures  are  being  developed  concurrently  with  the  compu¬ 
tational  effort. 

Three  different  approaches  to  the  solution  of  the  unsteady 
Euler  equations  were  examined  in  the  overall  program.  First, 
Eidelman  developed  a  two-dimensional  code  based  on  the 
Godunov  method  of  solution  [Ref.  1]  and  applied  the  code 
to  examine  unsteady  wave  propagation  in  ducts  [Ref.  2]  and 
the  process  of  port  opening  to  wave  rotor  passages  [Ref.  3],  A 
summary  of  Eidelman 's  work  is  given  in  Reference  4.  VJhile 
the  method  is  conservative  and  does  not  require  the  introduc¬ 
tion  of  artificial  viscosity,  the  extension  from  one  to  two 
dimensions  is  not  rigorous  when  shock  waves  are  present,  and 
computational  times  with  the  present  code  are  quite  long. 

The  extension  to  include  viscous  effects  would  require  a  separate 
treatment  of  the  boundary  layer. 

Second,  Mathur  developed  a  one-dimensional  code  based  on 
the  Random  Choice  Method  (RCM)  of  solution  [Ref.  5].  Somewhat 
similar  to  the  Godunov  method,  in  that  the  solution  is  based 
on  solving  the  Riemann  problem  within  each  grid  cell  at  each 
time  step,  the  RCM  approach  results  in  very  sharp  discontinui¬ 
ties  which  can  be  tracked  easily.  This  is  particularly  useful 
in  constructing  wave  rotor  cycles,  in  which  the  position  of 
the  gas-gas  interface  is  equally  as  important  as  the  position 


of  the  compression  and  expansion  waves.  The  code  is  there¬ 
fore  valuable  in  the  preliminary  design  process,  to  examine 
suitable  port  arrangements,  and  the  gas  properties  at  the 
ports,  for  a  given  task.  Unfortunately,  an  extension  to  two- 
dimensional  and/or  viscous  flow  can  not  be  made  rigorously. 

The  third  approach  was  followed  in  the  present  work.  The 
QAZID  method  for  compressible  inviscid  flow  computations 
developed  by  Verhoff  and  O'Neil  [Ref.  6]  was  implemented  to 
generate  a  one-dimensional  unsteady  Euler  code  with  the  goal 
of  evaluating  the  suitability  of  the  approach  for  wave  rotor 
applications.  Some  advantages  of  the  QAZID  method  were  recog 
nized  as  being  the  following: 

1.  The  method  is  based  on  the  use  of  characteristics. 

Such  methods  can  model  wave  propagation  accurately. 

2.  The  use  of  a  natural  streamline  coordinate  system 
eases  the  difficult  task  of  computing  with  two  and 
three  dimensional  grids. 

3.  The  equations  are  written  in  a  form  which  allows  a 
straightforward  extension  to  viscous  flows. 

4.  Codes  for  computing  internal,  steady,  two-dimensional 
flows  in  the  presence  of  shocks,  and  simple  internal 
viscous  flows,  have  been  generated  quickly  without 
significant  development  problems  [Ref.  6]. 

In  the  work  reported  herein,  a  one-dimensional  FORTRAN 
code  based  on  the  QAZID  method  was  developed  using  the  NFS 
IBM370-3033  computer  and  subsequently  exercised  on  the  shock 
tube  test  problem.  In  reporting  the  work,  first,  in  Section 
II  and  Appendix  A,  a  complete  account  is  given  of  the  deri¬ 
vation  and  non-dimensionalization  of  the  governing  equations. 


In  Section  III,  a  FORTRAN  code,  EULER-1,  is  described.  Its 
operation  on  the  NFS  computer,  and  the  listing  of  the  code, 
are  given  in  Appendix  B  and  Appendix  C,  respectively. 

Results  of  applying  the  code  to  the  shock  tube  test  problem 
are  given  in  Section  IV.  Difficulties  encountered  in  the 
implementation  of  the  method  and  additional  comments  are  given 
in  Section  V,  and  conclusions  are  given  in  Section  VT . 


II .  THE  QAZID  METHOD 


A.  OVERVIEW 

The  QAZID  method  uses  Riemann-like  variables  with  a  modi¬ 
fied  entropy  term  in  their  definitions,  to  express  the  Euler 
equations  in  a  natural  streamline  coordinate  system.  The 
resulting  partial  differential  equations  (PDE) ,  when  cast 
along  the  characteristic  trajectories  in  the  space-time 
domain,  reduce  to  a  system  of  ordinary  differential  equations 
(ODE)  which  may  be  solved  explicitly.  The  advantage  in  using 
the  modified  Riemann  variables  is  that  they  are  less  affected 
by  u; scontinuities  in  the  flow  than  are  the  standard  Riemann 
variables.  Since  the  equations  are  not  valid  across  discon¬ 
tinuities  which  cause  irreversible  losses,  such  discontinui¬ 
ties  must  be  located  and  treated  with  special  logic  in  the 
numerical  formulation. 

In  the  following  presentation,  the  reader's  familiarity 
with  the  method  of  characteristics  is  assumed. 

B.  DEVELOPMENT  OF  THE  EQUATIONS 

This  section  presents  the  governing  equations  and  outlines 
their  development.  A  rigorous  derivation  of  the  equations 
is  presented  in  Appendix  A. 

1 .  The  Coordinate  System 

The  natural  streamline  coordinate  system  fs,n,m),  is 
shown  in  Fig.  A-1  relati^'^e  to  a  fixed  rectangular  cartesian 


system  (x,y,z).  The  (s,n,m)  system  is  a  right-hand  orthogonal 
system  which  undergoes  curvilinear  translation  as  it  moves 
with  a  fluid  particle  along  a  streamline.  The  coordinate 
system  is  described  in  detail  in  Appendix  A. 

2 .  Variables 

The  modified  Riemann  variables,  or  "extended  Riemann 
variables"  [Ref.  6:p.  1],  are  defined  as 


Q  =  q  +  AS 
R  =  q  -  AS 


(1) 


where  q  is  the  velocity  magnitude,  A  is  the  speed  of  sound 
and  S  is  the  modified  entropy  defined  in  terms  of  pressure 
and  density  as 


S 


(2) 


The  modified  entropy  relation  is  the  result  of  defining  the 
entropy  change,  dS,  to  be  given  by 


dS  = 


(3) 


where  dQj^  is  the  heat  required  to  be  added  in  a  reversible 
process  between  the  same  end  states,  T  is  the  temperature  and 
r  is  the  ratio  of  specific  heats. 


3. 


Conservation  of  Mass 


By  applying  the  continuity  equation  to  a  differential 
stream  tube  of  variable  cross  sectional  area,  in  natural 
coordinates,  and  ignoring  all  third  order  anH  higher  deriva¬ 
tives,  one  obtains 


3  0 

3t 
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cos  ^ 
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(4} 


where  q  is  the  velocity,  o  is  the  density,  s  is  the  stream- 
wise  spatial  dimension,  and  6  and  $  are  the  flow  angles  as 
defined  in  Fig.  A-1  of  Appendix  A. 

4 .  Conservation  of  Momentum 

Applying  the  vector  form  of  the  momentum  conservation 
law  in  natural  coordinates,  the  equation  of  motion  for  inviscid 
flow  becomes 


'  r  9q  , 
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+  i^[p  qcos  0  If  +  P  q^'cos  0  ||  +  |1]  = 


5 .  The  Conservation  of  Energy 

In  the  absence  of  friction  and  heat  conduction,  and 
outside  of  discontinuities  (through  which  irreversible 
changes  consistent  with  mass,  momentum  and  energy  conservation 
are  permitted) ,  energy  conservation  is  equivalent  to  the 


statement  that  the  entropy  of  a  fluid  particle  does  not  change 
In  the  natural  coordinate  system,  therefore 


dS  ?S  ^  3S 

dt  =  ^  ^  ° 


Equation  (6)  will  not  be  valid  across  shock,  waves  or  contact 
surfaces  between  gases  having  different  states. 


6.  Transformation  to  a  Useful  Form 


Equations  (4)  and  (5)  are  first  expressed  in  terms 
of  the  primary  variables  q.  A,  S,  and  the  flow  angles.  Us 
the  definition  of  sound  speed  in  a  perfect  gas 


A  =  YP/P 


Eq.  (4)  becomes 


|A  +  q£A  +  i^A|f  +  q  +  cos  9 


=  0 


With  Eq.  (3),  the  equation  of  state  for  a  perfect  gas,  and 
the  first  law  of  thermodynamics,  Eq.  (5)  gives 
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Equations  (6),  (8)  and  (9)  constitute  the  system  of  P.D.E.'s 

which  describe  the  unsteady  isentroplc  flow  of  a  perfect  gas 
in  natural  coordinates.  Since  third  order  terms  were 
neglected  in  the  conservation  of  mass,  the  system  is  only 
accurate  to  second  order. 

The  system  of  P.D.E.'s  is  now  transformed  into  a  system 
of  O.D.E.'s  along  the  characteristic  trajectories.  In  general, 
if  w  is  a  function  of  (s,t),  then 

dw  =  ds  +  dt 

3s  3t 


or 


dw  _ 

dt  “  3s ^dt'  3t 


(10) 


where  (ds/dt)  describes  the  "characteristic"  direction  in  the 
(s,t)  plane  along  which  Eq.  (10)  gives  the  rate  of  change  of 
the  parameter  w.  After  the  appropriate  algebra  and  the  intro¬ 
duction  of  Eq.  (1),  the  final  system,  of  aquations  becomes 
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A  trajectory  in  the  space-time  domain  is  illustrated  in  F 


6t 


[6w  =  Wg-Wg  j 
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I  Aw  =  w  -w  ] 
Si  ^ 


Figure  1.  Solution  Procedure  in  the  Space-Time  Domain 


which  shows  an  infinitesimal  interval  of  space  between  two 
"nodes"  of  the  computational  mesh. 

For  the  change  in  w  from  A  to  B  along  the  trajectory 
with  slope  ' , 
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v'here  6w  denotes  the  change  due  to  time  at  a  fixed  location 
and  \w  denotes  the  change  due  to  displacement  at  a  fixed 
time,  for  the  characteristic  trajectory.  The  essence  of  the 
solution  procedure  is  to  calculate  the  change  in  the  varia¬ 
bles  at  each  spatial  node  during  the  differential  time 
interval,  so  that  the  step  to  the  next  time  level  can  be 
made.  In  other  words,  6w  must  be  calculated  at  each  node 
using 

_  _  t+5t  _ 

6w  =  -  Iw  +  /  Z  dt  (1 

t 

Since  all  the  information  at  time  level  t  is  known,  the 
position  of  A  can  be  determined  by  an  iterative  procedure 
based  on  a  and  Aw  can  be  calculated  by  interpolation.  The 
line  integral  can  be  evaluated  by  transforming  the  integral 
to  a  purely  spatial  integral  using  \  =  ds/dt.  Thus 

t+'t  B  s 

/  Z  dt  =  /  Y  ds  =  /  X 

t  A  '  s .  +  A  s 

1 

and  the  i.ntegral  can  be  evaluated  by  any  one  of  several 
numerical  methods.  Thus,  the  variables  at  each  node  can  be 
updated  in  time  using  Fq.  (14)  and  the  process  repeated. 


D.  DISCONTINUITIES  IN  THE  FLOW 

There  are  several  types  of  discontinuities  that  must  be 
considered.  They  are  gradient  discontinuities,  contact 
discontinuities  and  shock  waves. 

Gradient  discontinuities  are  characteristic  of  the  head 
and  tail  of  rarefraction  waves,  the  collision  of  two  shocks  and 
the  interaction  of  a  shock  and  a  contact  discontinuity. 

Across  such  a  gradient,  the  derivatives  of  velocity,  pressure 
and  sound  speed  are  discontinuous. 

A  contact  discontinuity  is  caused  by  the  interaction  of 
two  shocks  of  opposite  family  and  when  a  shock  overtakes  a 
shock  of  the  same  family.  Entropy  and  sound  speed  are  discon¬ 
tinuous  at  a  contact  discontinuity. 

Across  shock  waves,  velocity,  pressure,  density  and 
entropy  are  discontinuous  [Ref.  7]  . 

Since  Eq.  (11)  is  not  valid  across  discontinuities,  addi¬ 
tional  logic  is  required  in  the  numerical  procedure  to  model 
flows  which  contain  discontinuities.  The  method  presented  in 
Reference  6  for  making  a  correction  in  the  case  cf  shock  waves 
will  give  good  accuracy  in  problems  where  the  solution  is 
converging  to  a  steady  state  condition  but  will  not  give 
accurate  results  during  the  transient  portion  of  such  problems 
or  for  problems  with  only  unsteady  solutions.  In  the  case  of 
applications  to  the  wave  rotor,  it  would  certainly  be  necessary 
to  know  the  locations  of  the  contact  surfaces  between  the 


hot  and  cool  gases  as  well  as  the  locations  of  the  shock  and 


expansion  waves. 

The  methods  used  in  the  present  work  to  correct  for  dis¬ 
continuities  are  those  of  Moretti  [Ref.  7].  The  present 
discussion  will  be  limited  to  the  treatment  of  shock  waves. 
The  method  makes  use  of  the  analytical  relationship  between 
the  change  in  the  Riemann  variables  across  a  shock  and  the 
incoming  Mach  number  relative  to  the  shock  wave  (W)  which  is 
illustrated  in  Figs.  2  and  3.  The  relation  is  used  to 
determine  the  shock  speed  and  to  transform  the  problem  to 
a  steady  case  which  can  be  handled  using  normal  shock  rela¬ 
tions.  For  the  situation  depicted  in  Fig.  2  of  a  shock 
propagating  to  the  right  with  velocity  V  into  air  with 
velocity  q^,  and  with  the  high  pressure  side  to  the  left, 
where  A  denotes  the  left  side  of  the  shock  and  B  the  right, 


W 


B 
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If  the  pressure  and  density  are  non-di mensional ized  by  the 
values  on  the  low  pressure  side  of  the  shock,  the  change  in 
the  extended  Riemann  variable  Q  is  given  by 
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where 
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This  exact  relationship  can  be  approximated  by  the  polynomial 
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=  -2.7574  +  3.1573W  -  0.2863W 


(17 


as  illustrated  in  Fig.  4  over  a  shock  strength  range  of  1.0 
to  4 . 0 .  It  should  be  noted  that  if  the  high  pressure  side 
were  to  the  right,  as  in  Fig.  3,  Eq.  (15)  would  become 


Vv 


“a 


-  V 

A  s 


and  the  left  side  of  Eq.  (16)  would  become 
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The  procedure  then,  is  to  measure  the  change  in  Q  across  an 
interval  where  a  shock  is  known  to  exist.  This  value,  call 
it  is  used  in  Eq.  (17)  to  get  an  approximation  for  the 

corresponding  value  of  W  and  then  the  exact  value  of  AQ,  call 
it  AQ^,  is  calculated  using  Eq.  (16)  .  If  the  exact,  value  of 
AQ  Is  not  equal  to  the  value  measured  across  the  shock,  a 


new  value,  AQ,  is  calculated  according  to 


and  entered  into  Eq.  (17)  to  obtain  another  value  of  W.  When 
the  exact  value  of  AQ  calculated  from  Eq.  (16)  equals  the 
measured  value  of  AQ,  W  is  )cnown  to  be  correct  and  the  normal 
shock  relations  can  be  used  to  calculate  the  values  at  node  A 
Also,  Vg  can  be  calculated  from  Eq .  (15)  and  used  to  track 
the  shock  during  the  next  time  interval. 


III. 


DESCRIPTION  OF  FORTRAN  PROGRAM  " EULER- 1 


A.  machine  and  language 

The  EULER-1  program  is  written  in  VS  FORTRAN  and  runs  on 
an  IBM  3033,  System  370  computer,  however  the  code  is  simple 
and  small  enough  to  enter  and  run  on  a  mini  or  micro  cnmpiiter 
in  a  Basic  language  if  one  is  willing  to  accept  significantly 
longer  run  times.  Table  1  at  the  end  of  this  section  is  a 
summary  of  the  editable  parameters  and  their  effect  on  the 
program. 

B.  CONVENTIONS  AND  BASIC  STRUCTURE 

EULER-1  is  a  first-order  one-dimensional  code  using  an 
evenly  spaced  numerical  grid.  All  values  are  double  precision 
except  those  used  in  the  graphics  routines  which  are  rounded 
to  single  precision. 

Each  subroutine  has  its  own  variables  space.  Jn  other 
words,  variables  are  not  shared  in  common  throughout  the 
program  but  must  be  passed  to  the  subroutine  being  called 
by  the  calling  routine.  In  all  cases,  however,  the  variables 
have  the  same  name  in  both  the  called  and  the  calling  routine 
so  there  should  be  no  confusion.  This  was  done  so  that 
arrays  could  be  dimensioned  at  execution  time. 

The  program,  depicted  in  Fig.  5,  is  structured  around  a 
main  routine  which  serves  as  a  user  input  area,  sets  up  the 
problem  and  calls  five  subroutines  to  solve  the  problem  and 


call 

"EXACT" 


output  the  solution.  The  five  subroutines  are  "TIME,"  "TRAK," 
"SWEEP,"  "JUMP"  and  then  one  of  several  output  routines  depend¬ 
ing  on  user  desires.  Output  options  include  two  types  of 
graphical  displays,  "PLOT"  and  "EXACT,"  and  one  tabular 
listing  routine,  "LIST."  When  "PLOT"  is  selected,  "BORDER" 
is  automatically  called  to  set  up  the  plotting  area.  There  are 
two  other  subroutines,  "RSHOCK"  and  "LSHOCK,"  which  are  called 
by  "SWEEP"  as  needed. 

C.  SUBROUTINE  DESCRIPTIONS 

In  general,  each  subroutine  begins  with  a  heading  followed 
by  variable  definitions  v;here  appropriate,  variable  declara¬ 
tions  and  array  dimensioning.  Most  variables  are  defined  in 
the  "MAIN"  routine  and  only  those  variables  which  were  not 
are  defined  in  the  subroutines  where  they  are  used. 

1.  The  "MAIN"  Routine 

This  routine  forms  the  main  structure  of  the  program 
and  includes  a  heading,  an  extensive  list  of  variable  defini¬ 
tions,  a  user  input  area,  the  necessary  statements  to  make 
the  initial  value  assignments  to  variables,  and  the  call 
statements  for  the  various  subroutines. 

The  input  area  is  those  lines  of  the  program  (145-175) 
where  the  user  edits  the  program  to  establish  the  initial 
conditions,  mesh  size,  termination  criteria  and  output  options 

The  initial  conditions  which  m.ay  he  modified  are  the 
pressure,  temperature  and  density  ratios  across  the  diaphragm, 
the  initial  velocity  of  the  fluid  on  each  side  of  the 


diaphragm,  and  the  value  of  gamma  (which  must  be  the  same  for 
both  sides) . 

All  velocities  are  non-dimensionalized  by  the  initial 
sound  speed  on  the  low  pressure  side  of  the  diaphragm  and 
pressures  and  densities  are  non-dimensionalized  by  their 
initial  values  on  the  low  pressure  side  of  the  diaphragm. 

The  mesh  size  may  be  set  at  any  odd  number  and  the 
arrays  in  the  user  input  area  must  be  dimensioned  as  such. 

The  criteria  for  program  termination  is  the  number  of 
time  steps  computed,  which  the  user  selects. 

The  output  options  are  determined  by  the  value  of  the 
variable  GRAPHS.  A  value  of  zero  results  in  a  call  to  "LIST" 
which  writes  to  file  9  on  the  -.user’s  permanent  disk,  a  tabular 
listing  of  the  variable  arrays  and  discontinuity  locations. 

A  value  of  1  results  in  a  call  to  "BORDER"  and  "PLOT"  which 
create  a  plot  of  the  pressure,  density,  velocity  and  entropy 
distributions  as  in  Fig.  7.  A  value  of  2  results  in  a  call 
to  "EXACT"  which  creates  a  plot  of  the  density  distribution 
compared  to  the  exact  solution  as  in  Fig.  8.  The  exact  values 
to  be  plotted  must  be  entered  in  the  user  input  area,  they  are 
not  calculated  by  the  program.  A  more  detailed  description 
of  the  various  outputs  is  found  under  the  appropriate  subrou¬ 
tine  description.  The  frequency  with  which  output  is  created 
is  controlled  by  the  variable  SKIP  which  is  the  number  of  time 
steps  between  calls  to  output  routines. 


The  "TIME"  Routine 


The  maximum  allowable  time  step  is  governed  by  the 
CFL  condition.  Simply  stated,  this  means  that  the  time  step 
must  be  small  enough  so  that  the  characteristic  trajectories 
remain  within  one  spatial  interval  during  the  time  interval. 
The  minimum  time  step  is  calculated  by  computing 

Delt  =  H/ABS[Q+A] 

at  every  node,  where  H  is  the  spatial  interval,  and  selecting 
the  minimum  value  of  Delt. 

3 .  The  "TRAK"  Routine 

Shock  locations  at  the  unknown  time  level  are  deter¬ 
mined  by  computing  the  shock  speed  at  the  known  time  level, 
as  outlined  in  Section  2,  and  multiplying  by  the  time  inter¬ 
val  computed  by  "TIME."  The  time  step  is  reduced,  if  neces¬ 
sary,  to  limit  shock  travel  to  one  spatial  interval.  The  node 
imm.ediately  to  the  right  of  the  shock  (upstream)  is  flagged 
for  later  use  by  "SWEEP"  and  "JUMP." 

In  calculating  the  shock  speed,  the  assumption  is  m.ade 


that  the  conditions  immediately  adjacent  to  the  shock  are  the 
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time  level  for  the  values  of  Aw,  calculation  of  z,  and 
integration  of  z  for  each  node. 

To  calculate  Aw,  the  assumption  is  made  that  the 
characteristic  trajectories  are  straight  lines.  An  initial 
guess  of  the  characteristic  slope  (A)  is  made  and  by  forcing 
the  characteristic  to  pass  through  the  node  at  point  B  in 
Fig,  1,  As  =  A  xdelt.  The  values  of  q  and  A  at  point  A  are 
found  by  interpolation  and  the  slope  of  the  characteristic 
is  then  calculated  using  the  values  of  q  and  A  at  point  A. 
This  slope  is  compared  with  the  initial  guess  and  if  the  two 
are  not  in  agreement,  the  procedure  is  iterated  until  they 
converge.  Once  As  is  accurately  determined,  the  values  of  Aw 
can  be  calculated  by  interpolation.  Linear  interpolation  is 
used  in  the  present  version  of  EULER-1.  This  procedure  is 
carried  out  for  each  characteristic  at  each  node. 


i-1  i  i+1 


Figure  6.  "F.SHOCK"  Case 

A  deviation  from  the  above  procedure  is  necessary  in 
the  case  where  a  shock  exists  in  the  spatial  interval  to  the 
left  or  right  of  the  current  node.  When  a  shock  exists  to 


the  right,  and  the  flow  is  subsonic,  as  in  Fig.  6,  care  must 
be  taken  not  to  interpolate  for  point  3  based  on  the  slope 
between  points  T+1  and  T,  which  would  not  be  accurate  due 
to  the  discontinuity.  In  such  a  case  "RSHOCK"  is  called  and 
the  interpolation  is  based  on  the  slope  between  point  T  and 
I-l.  Similarly,  when  the  shock  is  to  the  left  of  the  current 
node,  "LSHOCK"  is  called  and  the  interpolation  is  based  on 
the  slope  between  I+l  and  I. 

The  calculation  of  z  includes  the  calculation  of 
spatial  derivatives  of  q  and  A  for  characteristics  1  and  3 
of  Fig.  6.  In  general,  the  derivatives  associated  with 
characteristic  3  are  forward  differenced  and  those  associated 
with  characteristic  1  are  backward  differenced  in  keeping 
with  the  principle  of  domain  of  dependence.  If  the  flow  is 
supersonic  or  if  a  shock  exists  to  the  right  of  the  current 
node,  all  derivatives  are  backward  differenced.  If  a  shock 
e.xists  to  the  left  of  the  current  node,  all  derivatives  are 
forward  differenced.  Once  the  derivatives  are  known,  z  is 
calculated  from  Eq .  (11)  using  the  average  value  of  A  between 

points  1  and  I  or  3  and  I  as  appropriate.  Since  the  deriva¬ 
tives  are  linear,  this  results  in  an  average  value  of  z 
over  the  same  interval. 

The  integration  of  z  is  transformed  from  a  time  to  a 
spatial  integration  as  described  in  Section  IT  and  the 
trapezoidal  rule  is  iised  to  carry  out  the  integration.  In 
EULER-1,  this  has  been  done  in  one  step  using  the  average 
''•alue  of  z  described  above. 


Equation  (14)  is  then  solved  by  addition  of  aw  and 
the  results  of  the  integration,  and  the  resulting  6w  i s 
stored  until  all  nodes  are  calculated.  After  all  nodes  have 
been  calculated,  the  variable  arrays  (w)  are  updated  for  the 
next  time  interval. 

5 .  The  "LSHOCK"  Routine 

As  discussed  above,  the  "LSKOCK"  routine  modifies 
the  basic  EULER-1  procedure  at  an  interior  node  when  a  shock 
exists  to  the  left  of  the  node.  Interpolation  of  quantities 
to  the  left  of  the  node  are  computed  based  on  the  derivatives 
of  the  quantities  to  the  right  of  the  node.  Although  this 
assumes  that  the  derivatives  do  not  change  between  adjacent 
spatial  intervals,  it  is  necessary  to  avoid  taking  derivatives 
across  discontinuities  in  the  flow. 

6 .  The  "RSHOCK"  Routine 

Similar  to  "LSHOCK,"  "RSHOCK"  bases  interpolations 
of  quantities  to  the  right  of  the  node  on  the  derivatives  of 
the  quantities  to  the  left  of  the  node  when  a  shock  exists 
to  the  right  of  the  node. 

7 .  The  "JUMP"  Routine 

The  "JUMP"  routine  is  used  to  calculate  the  conditions 
downstream  of  a  stationary  shock  as  described  in  Section  II. 

If  a  shock  is  known  to  have  crossed  a  spatial  node  during  a 
time  interval,  which  is  known  once  "TRAK"  has  been  called  for 
that  time  interval,  the  entire  "SWEEP"  sequence  is  skipped  for 
the  node  which  was  crossed  by  the  shock  and  the  conditions  at 


that  node  are  determined  using  the  normal  shock  relations.  No 
that  the  node  in  q\iestion  is  the  node  downstream  of  the  sta¬ 


tionary  shock.  As  in  "TRACK,"  it  is  assumed  that  the  condi¬ 
tions  at  the  nodes  upstream  and  downstream  of  the  shock  are 
the  same  as  the  conditions  immediately  adjacent  to  the  shock. 

8 .  The  Output  Routines 

There  are  four  subroutines  which  produce  various  types 
of  output  as  the  user  desires.  "BORDER"  and  "PLOT"  produce 
a  graphical  presentation  of  the  pressure,  density,  velocity 
and  entropy  distributions  at  selected  time  levels  as  seen 
for  example,  in  Fig.  7.  "EXACT"  produces  a  plot  of  the  den¬ 
sity  distribution  at  one  selected  time  interval  and  compares 
it  with  the  exact  solution  as  shown  in  Fig.  8.  At  selected 
time  levels,  "LIST"  produces  a  tabular  listing  of  the  Riemann 
variables,  modified  entropy,  pressure,  density  and  velocity 
distributions,  elapsed  time,  shock  speed  and  discontinuity 
locations.  The  listing  is  written  to  the  user's  permanent 
storage  disk. 

When  the  value  of  GRAPHS  is  set  equal  to  one  in  the 
"MAIN"  routine,  "BORDER"  is  called  once  to  set  up  the  plot 
axis,  labels,  and  headings.  "PLOT"  is  called  every  SKIP 
time  steps  to  draw  the  four  distribution  curves. 

When  the  value  of  GRAPHS  is  set  equal  to  two,  "EXACT" 
is  called  every  SKIP  time  steps  to  plot  the  density  distri¬ 
bution  compared  with  the  exact  solution  computed  at  six 
points  of  interest.  The  points  are  the  two  end  points,  which 


->  O 


are  simply  the  initial  conditions,  the  head  and  tai],  of  the 
rarefraction  wave,  the  point  just  left  of  the  contact  surface 
and  the  point  just  left  of  the  shock.  The  spatial  locations 
of  these  points  are  computed  by  "EXACT"  based  on  the  elapsed 
time  and  the  known  values  of  the  wave  velocities  entered  in 
the  "MAIN"  routine.  The  exact  values  of  the  densities  at  these 
points  must  also  be  entered  in  the  "MAIN"  routine. 

When  the  value  of  GRAPHS  is  set  equal  to  zero,  the 
tabular  listing  as  described  above  is  sent  to  the  user's 
disk.  No  graphical  output  is  created. 

D.  CAPABILITIES  AND  LIMITATIONS 

The  present  version  of  EULER-1  is  set  up  to  solve  a  shock 
tube  problem  with  a  single  centered  diaphragm.  Boundary 
conditions  for  the  ends  of  the  tube  have  not  been  incorporated 
so  the  problem  must  be  stopped  before  the  waves  reach  the  end 
of  the  tube. 

The  left  side  of  the  diaphragm  is  the  high  pressure  side. 
The  program  will  not  run  with  the  right  side  as  the  high 
pressure  side  without  changes  to  some  of  the  shock  correction 
logic. 

The  user  may  select  any  odd  number  of  grid  points  limited 
only  by  the  amount  of  memory  available. 

The  program  tracks  shock  v;aves  and  makes  shock  jump  calcu¬ 
lations  at  the  appropriate  locations  but  the  program  can  also 
run  without  the  shock  tracking  feature  and  jump  calculations 
if  so  desired  with  some  smearing  of  the  shock  discontinuity 
and  complete  loss  of  entropy  change  across  the  shock. 
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TABLE  1 


LIST  OF  EDITABLE  PARAMETERS 


Parameter 


Function 


N  Number  of  grid  points 

M  Controls  which  discontinuities  are 

tracked 

1  =  Shocks  only 

2  =  Shocks  and  contact  surfaces 


GRAPHS  Controls  the  form  of  the  output 

0  =  Tabular  listing 

1  =  Pressure,  density,  velocity, 

entropy  plot 

2  =  Exact  solution  comparison  for 

density 

SKIP  Number  of  time  steps  between  output 

calls 


JSTOP 

TRI 

PRI 

DRI 

QLI 

QRI 

G 


Number  of  time  steps  calculated 
Initial  temperature  ratio 
Initial  pressure  ratio 
Initial  density  ratio 

Initial  velocity  left  of  the  diaphragm 
Initial  velocity  right  of  the  diaphragm 
Ratio  of  specific  heats 


EULER-1  was  tested  on  the  Riemann  shock  tube  problem. 

The  initial  conditions  and  the  courseness  of  the  computa¬ 
tional  mesh  were  varied.  Additionally,  one  run  was  made 
without  the  shock  tracking  and  correction  features. 

Figures  7  and  8  show  the  results  for  initial  pressure 
and  density  ratios  of  5,  uniform  temperature  and  with  the  air 
initially  at  rest,  using  a  grid  of  101  points.  The  exact, 
analytically  predicted  conditions  are  also  shown  in  Fig.  8. 

It  can  be  seen  that  all  wave  velocities  are  correctly  computed 
and  the  shock  is  defined  within  one  interval .  The  rarefrac- 
tion  waves  are  slightly  smeared  and  the  contact  surface  is 
greatly  smeared.  A  slight  transient  instability  to  the  right 
of  the  diaphragm  is  noticeable,  particularly  in  the  plots 
of  pressure  and  velocity. 

Figures  9  and  10  show  results  for  initial  pressure  and 
density  ratios  of  1.3.  The  shock  is  still  well  defined  in 
the  correct  interval  and  the  contact  surface  is  not  as  badly 
smeared  as  in  Fig.  8.  There  is  a  loss  of  accuracy  however, 
with  respect  to  the  tail  of  the  rarefraction  wave. 

Figure  11  shows  the  results  of  the  original  problem  with 
a  grid  of  51  points.  As  might  be  expected  with  a  coarser 
mesh,  the  discontinuities  are  less  sharply  defined  although 
it  is  clear  that  the  wave  velocities  are  accurately  computed. 


Figure  12  shows  the  results  for  the  original  problem  with 
a  grid  of  901  points.  It  can  be  seen  that  the  EULER-1 
solution  closely  approximates  the  exact  solution. 

Figures  13  and  14  show  the  result  of  running  the  original 
problem  without  the  shock  tracking  and  correcting  features. 
Note  that  the  shock  position  is  incorrectly  computed  and  it 
is  smeared  slightly.  Also  note  in  Fig.  13,  the  lack  of  any 
change  in  modified  entropy  across  the  shock. 

Run  time  for  the  101  point  mesh  is  0.00049  seconds  per 
node  per  time  step.  The  results  for  Fig.  7  took  50  time 
steps  for  a  total  time  of  2.46  seconds.  For  the  901  point 
mesh,  the  run  time  decreased  to  0.00023  seconds  per  node  per 
time  step  due  to  the  less  frequent  use  of  the  shock  tracking 
and  correction  routines  per  node  calculated.  The  results  for 
Fig.  12  took  470  time  steps  for  a  total  time  of  71.85  seconds 
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Figure  7. 


EULER-1  Results  for  Pressure  Ratio  5  with 
Shock  Tracking 
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Figure  8.  Exact  Solut 
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Figure  11..  Course  Mesh  Solution  for  Pressure  Ratio 
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Figure  12.  Fine  Mesh  Solution  for  Pressure  Ratio 
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Figure  13. 

EULEP-1  Results  for  Pressure 
without  Shock  Tracking 
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Figure  14,  Exact  Solution  Comp 
without  Shock  Track 


V.  DISCUSSION 


A.  RESULTS 

The  transient  instability  noted  in  the  results  was  not 
investigated  since  the  amplitude  was  small  and  the  effect 
would  be  of  no  consequence  in  most  applications. 

At  low  pressure  ratios,  the  velocity  of  propagation  of 
the  tail  of  the  rarefraction  wave  approached  sound  speed 
and,  in  these  cases,  the  EULER-1  code  did  not  accurately 
compute  its  location.  The  reason  for  this  was  not  investi¬ 
gated  here;  however,  due  to  the  low  pressure  ratios  at  which 
wave  rotors  can  operate  /  this  is  of  interest  and  should  be 
investigated  in  future  work. 

B.  SPECIAL  CONSIDERATIONS 

1 .  Modified  E.ntropy 

In  the  QAZID  method,  and  subsequently,  in  the  EULER-1 
code,  the  variable  S,  which  here  has  been  referred  to  as  the 
"modified  entropy"  and  which  is  referred  to  in  Reference  6  as 
simply  the  "entropy,"  does  not  behave  thermodynamically  as 
one  would  expect  entropy  to  behave.  As  a  fl\iid  crosses  a 
shock  wave,  the  modified  entropy  of  the  fluid  decreases. 

This  is  the  result  of  the  definition  of  Eq.  (3) ,  repeated  here 


The  use  of  the  word  "entropy"  and  symbol  "S"  for  this  varia¬ 
ble  is  confusing  since  the  sign  of  the  change  in  "S"  is  in 

direct  conflict  with  well-established  conventions  in  thermo- 

nd 

dynamics.  For  example.  Corollary  6  to  the  2  Law  of 
Thermodynamics  [Ref.  8;p.  88]  would  have  to  be  changed  to 
read,  "The  'modified  entropy*  of  an  isolated  system  decreases 
or  in  the  limit  remains  constant."  At  a  minimum  perhaps,  the 
symbol  §  to  denote  the  "modified  (negative  and  scaled)  entropy" 
would  be  helpful.  In  practice  of  course,  the  change  in  the 
conventionally  defined  entropy  can  be  recovered  from  the  change 
in  the  modified  entropy,  simply  by  multiplying  by  (-y) • 

2 .  Moretti's  Methods 

Incorporating  Moretti's  methods  of  handling  discon¬ 
tinuities  [Ref.  7]  into  the  EULER-1  code  required  special 
care . 

There  is  a  fundamental  difference  in  procedure  in 
that  the  EULER-T  code  is  based  on  the  high  pressure  side 
being  on  the  left  whereas  Moretti's  formulations  are  all 
based  on  the  high  pressure  side  being  on  the  right.  Although 
essentially  a  matter  of  bookkeeping,  it  is  an  area  of  poten¬ 
tial  confusion. 

A.nother  difference  is  that  the  Riemann  variables  used 
by  Morreti  are  not  the  Riemann  variables  used  in  the  QAZlD 
method.  This  can  lead  to  difficulty.  In  fact,  in  the 
application  of  Moretti's  shock  tracking  scheme  to  the  EULER- 1 
code,  the  pressures  and  densities  have  to  be  non-dimensional ized 


by  their  values  on  the  low  pressure  side  of  the  shock  in  order 
to  obtain  Eq.  (14),  and  for  the  procedure  to  work. 


C.  SUITABILITY  FOR  WA\^  ROTOR  APPLICATIONS 

Further  work  is  necessary  before  the  suitability  of  the 
QAZID  method  for  wave  rotor  applications  can  be  fully  evalu¬ 
ated.  The  following  areas  need  to  be  addressed. 

1 .  Boundary  Conditions 

In  the  EULER-1  code,  boundary  conditions  at  the  ends 
of  the  passage  have  not  been  incorporated.  The  code  calculates 
the  changes  of  the  interior  nodes  only,  skipping  the  two  end 
nodes.  In  particular,  solid  wall  boundary  conditions  and 
open-end  boundary  conditions  must  be  incorporated.  No  particu¬ 
lar  problems  are  expected  in  the  boundary  conditions  themselves 
However,  additional  logic  may  be  necessary  in  the  handling  of 
discontinuities,  such  as  the  reversal  of  conditions  to  the 
left  and  right  of  a  discontinuity  after  reflection  from  a 
boundary.  The  additional  logic,  which  may  be  considerable, 
is  warranted  by  the  simplicity  of  the  QAZlD  method. 

2 .  Contact  Discontinuity 

No  special  attention  is  given  to  the  contact  discon¬ 
tinuity  in  the  EULER-1  code  and  hence  the  discontinuity  is 
smeared.  In  a  wave  rotor  application,  this  discontinuity 
will  have  to  be  sharp.  Moretti's  methods  again  appear  to  be 
applicable  here  and  the  additional  logic  needs  only  to  be 
formulated  and  incorporated  into  the  EULER-1  code. 


VI.  CONCLUSIONS 


The  development  of  the  EULER  equations  in  the  natural 
streamline  coordinates  using  extended  Riemann  variables  was 
reviewed  in  detail.  The  QAZID  solution  method,  with  the 
addition  of  shock  tracking  features,  was  implemented  in  a 
first  order,  one-dimensional  FORTRAN  code  with  the  intent  of 
evaluating  the  method's  suitability  for  wave  rotor  applica¬ 
tions.  The  code  (EULER-1)  was  tested  on  the  shock  tube 
problem  with  good  results.  The  incorporation  of  boundary 
conditions,  an  improvement  in  the  contact  discontinuity 
definition  and  the  addition  of  an  area  variation  term  for 
quasi-one-dimensional  modeling  are  considered  to  be  neces¬ 
sary  before  the  QAZID  method  can  be  accepted  as  being  suitable 
for  wave  rotor  applications.  The  additional  logic  required 
for  these  extensions  may  be  considerable  but  the  development 
is  warranted  by  the  overall  simplicity  of  the  QAZID  method, 
and  its  straightforward  extension  to  viscous,  multi-dimensional 


APPENDIX  A 


DERIVATION  OF  THE  GOVERNING  EQUATIONS  FOR 
THREE-DIMENSIONAL  INVISCID  FLOW  (WITH  AREA  CHANGE) 

The  governing  equations  of  a  3-D,  inviscid,  compressible 
unsteady  flow  are  derived  here  in  a  natural  streamline 
coordinate  system.  The  equations  are  then  recast  along 
characteristic  trajectories  in  the  (s,t)  plane  and  expressed 
in  extended  Riemann  variables,  reducing  the  system  of  par¬ 
tial  differential  equations  to  a  system  of  ordinary  differ¬ 
ential  equations. 

A.  DEFINITION  OF  VARIABLES 

A  Cross-sectional  area  of  a  differential  stream 

tube 

A  Speed  of  sound 

Cp  Specific  heat  at  constant  pressure 

Specific  heat  at  constant  volume 

da  Normal  vector  for  differential  area  of  control 

surface 

e  Specific  internal  energy 

i  Unit  vectors  in  the  s.  n,  and  m  directions 

s ,  n,m 

P  Static  pressure 

Q  Modified  Riemann  variable 

Reversible  heat  transferred 

q  Velocity  magnitude 

r  Position  vector 


R 


Modified  Riemann  variable 


R|2  Gas  constant 

S  Modified  entropy 

T  Static  temperature 

t  Time 

u  Velocity  relative  to  a  standing  shock  wave 

V  Specific  volume 

vol  Differential  element  of  volume 

V  Velocity  vector 

W  Incoming  Mach  number  relative  to  a  standing 

shock  wave 

p  Density 

V  Ratio  of  specific  heats 

9,'+!  Flow  angles  with  respect  to  reference  coordinate 
planes 

V  Vector  operator 
B.  THE  COORDINATE  SYSTEM 

The  natural  streamline  coordinate  system  (s,n,m)  is 
shown  with  respect  to  a  fixed  rectangular  cartesian  system 
(x,y,z)  in  Fig,  A1 .  The  system  is  a  curvilinearly  translating, 
right  hand  orthogonal  system  which  translates  with  a  fluid 
particle  along  a  streamline  such  that  the  's'  coordinate  is 
measured  in  the  direction  of  the  flow.  The  'n'  coordinate 
direction  always  lies  in  the  plane  defined  by  the  's'  coor¬ 
dinate  direction  and  the  fixed  'y'  axis.  The  'm'  direction  is 
norma]  to  the  {s,n)  coordinate  surface.  The  flow  angles,. 

9  and  i ,  are  defined  as  shown  in  Fig.  A1 . 


Figure  A1 .  The  Natural  Coordinate  System 

C.  VECTOR  OPERATORS  IN  THE  NATURAL  COORDINATE  SYSTEM 
1.  V(  ) 

If  dr  is  the  position  vector  of  a  fluid  particle 
d?  •  V(  )  =  d(  ) 

In  natural  coordinates,  Eq.  (Al)  becomes 


[isds  +  i^dn  +  i^dm]-[(  )is  +  <  <  ^^m^ 


Us  iLUln  ^  M-Um 


8s 


9n 


9m 


By  inspection. 


2. 


V7(  ) 


Since 

V  =  q  is 


using  Eq.  (A2) 
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3.  7-V 

From  Eq.  (A2)  with  V  =  qX^ 


•  V  = 


v-v 


3s 


(A4 


4.  (V-V)V 


From  Eq.  (A3)  with  V  =  qi^ 
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(vv)v  =  q  ^  =  q 


3(qis) 


3i 


3s 


r  3q? 

~  3s^s  3s 


(A5 


The  change  in  the  unit  vector  i^  with  respect  to  s  is 
derived  below  and  is  illustrated  in  Figs.  A2  through  A6.  (A 
three-dimensional  model  is  very  helpful  in  visualizing  the 
vector  geometry.) 

Figure  A2  illustrates  the  natural  coordinate  system 
translating  along  a  streamline  from  point  A  to  point  B  in 


E 


the  inertial  cartesian  coordinate  system.  Figure  A3  is  the 

superposition  of  the  i^  unit  vector  from  points  A  and  B  of 

Fig.  A2 .  The  intent  here  is  to  express  the  change  in  the 

unit  vector  i  in  terms  of  components  along  the  original 

(s,n,m)  directions  shown  at  point  A  in  Fig.  A3.  Vector  OA 

in  Fig.  A3  is  the  original  i  direction  and  vector  OB  is 
di  ® 

g 

i_  +  ds.  Points  C  and  E  are  the  projection  of  points  B 

3  o  S 

and  A  in  the  (x,z)  plane  respectively.  Point  D  is  the  projec¬ 


tion  of  point  C  onto  the  line  OE. 

Figure  A4  is  the  plane  formed  by  the  OB  vector  and  the 
Y  axis.  The  length  of  OB  is  unity  by  definition.  The  length 


of  BC  is  sin(9+d0) 


sin(fi+d0)  =  sin  9  cos  dO  +  cos  6  sin  d0 


which,  since  dB  is  a  small  angle,  is 


sin(0+d9)  =  sin  9  +  cos  0  d9 


The  length  of  OC  is  cos(0+d6] 


cos(0+d0)  =  cos  0  cos  d9  -  sin  0  sin  d0 


which,  since  d0  is  a  small  angle,  is 


cos(9+dft)  =  cos  0  -  sin  0  dB 


Figure  A5  is  the  (x,z)  plane  in  which  the  angle  $  is 


defined.  The  i  direction  is  shown  at  point  C.  The  length 
m  ^ 

of  CD  is  OC  sindg  which,  with  Eg.  (A7)  is 

[cos  0  -  sin  0  d9]sin  d^) 

which,  since  dcf  is  a  small  angie,  is 

cos  9  dp  -  sin  0  d0  dp 

Dropping  the  second  order  term 

CD  =  cos  0  dp  (A8) 

The  length  of  OD  is  OC  cos  dp  which,  since  dp  is  a  small  angle, 
i  s  OC. 

Figure  A6  is  the  plane  formed  by  the  OA  vector  and 
the  y  axis.  GD  is  the  projection  of  BC  and  since  both  are 
vertical  lines,  BC  and  GD  have  the  same  length,  which  has 

^  /N 

been  determined  in  Eg,  (A6) .  The  i  and  i  directions  are 

s  n 

shown  at  point  A.  The  lengths  to  be  determined  are  AF  and  GF, 
since  these  are  the  components  of  AB  in  the  s  and  n  directions 
For  small  angles,  OD  was  shown  to  be  equal  to  OC  and  GD  =  BC , 
so  that  the  angle  GOD  is  equal  to  angle  BOC  =  0+d9.  Angle  GOF 
is  therefore  d9,  the  length  of  OG  is  unity  and  the  length 
of  GF  is  sin  dO,  which  for  small  angles  i.s  d0 ,  and  it  is  in 
the  i  direction.  The  length  of  AF  is  OA-OF.  OA  is  unity 


by  definition  and  OF  is  cos  d8,  therefore  AF,  the  component  of 


AB  in  the  i^  direction  is  equal  to  1  -cos  d6,  which,  since 
d0  is  a  small  angle,  is  zero. 

The  component  of  AB  in  the  i^  direction  is  GF  or 

de  i^  (A9) 


The  remaining  component  of  AB  in  the  i  direction  is 
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CD  or 
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Therefore,  the  vector  AB  can  be  written  as 
3i  .  . 


3s 


ds  =  de  +  cos  0  dd)  i 
m  ^ 
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ds  i  +  cos  0  ds  i 
dS  n  Ds  m 


and,  dividing  by  ds,  finally 
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Note  that  Fig.  A3  can  be  considered  also  to  describe 

/\ 

a  change  in  the  unit  vector  i^  at  a  fixed  location,  with 

respect  to  time.  In  this  case,  the  vector  AB  is  equal  to 
Di 


with  Eq.  (All),  Eq.  (A5)  becomes 
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D.  CONSERVATION  OF  MASS 


In  the  quasi-one-dimensional  differential  stream  tube, 
shown  in  Figure  A7 ,  where  is  always  in  the  direction  of 
V,  0,  q,  and  the  cross  sectional  area  A  are  known  at  the  center 
of  the  element  and  p  and  q  are  assumed  constant  on  any  given 
cross  section. 
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Figure  A7 .  Differential  Stream  Tube  of  Variable  Cross 
Section 


The  statement  of  continuity  is 


The  change  in  th<=‘  mass 
within  the  control  volume 
with  respect  to  time 


The  net  influx  of 
mass  across  the 
control  surface 


since  the  volume  is  constant,  the  left  side  becomes 


A  as  If 


For  the  stream  tube,  mass  enters  and  leaves  only  from  the 
ends  so  that  the  right  hand  side  can  be  written  as 


,  3p  ds,  r  3q  ds,  9A  ds,  ,  3p  ds,  .  dSi  , 

-a?  — >  '‘J-st  — >  ‘*-3?  T'  '  !'>  -^3?  — >  >'3  -^3?  + 


which  on  expansion  becomes 


-  pq  II  ds  -  qA  If  ds  -  pA  ||  ds 


3p  9q  3A  ,  3 

^  ^  I? 


On  dropping  the  higher  order  terms  and  combining  with 
the  left  hand  side, 


3p  ,  1  3a  ,  3p  .  9q  n 


3t 


3s 


3s 


3A  ds 
3s  2 


(A14) 


Equation  (A14)  is  the  form  of  the  continuity  equation 

which  is  required  to  model  a  flow  as  being  one-dimensional 

1  3A 

with  area  change.  In  general,  however,  —  must  be  expressed 
in  terms  of  q,  0  and  4> . 

With  reference  to  Fig.  A8 ,  the  change  in  the  cross- 
sectional  area  of  the  stream  tube  can  be  written  as 


3A 


=  [dm  + 


3  dm 


ds) [dn  + 


ds  1  - 


ds 


3dn 


A 


(A15 


+ 


The  net  force 
on  the  con¬ 
trol  surface 


The  change  in  the 
roanentum  of  the  fluid 
in  the  control  volume 
with  respect  to  time 


The  net  influx 
of  momentim 
across  the  con¬ 
trol  surface 


J 


+ 


The  net 
body 
force 

on  the  fluid 


Assuming  the  fluid  is  inviscid,  the  only  forces  acting  on  the 
control  surface  are  pressure  forces.  In  the  absence  of  any 
body  forces  (gravity,  electromagnetic,  etc.),  the  statement 
becomes 


// /  dvol  =  -  //  VpV'da  -  //  P  da  =  0  (A2 

With  Gauss'  Theorem,  Eq.  (A25)  becomes 
///  ■|^[pV]dvol  +  /// V  •  [VpV]  dvol  +  ///  VP  dvol  =  0 


or 


///{|^[pV] 


+  V’[VpV]  +  VPldvol 


0 


(A2' 


Since  Eq.  (A26)  is  true  for  any  volume,  the  integrand  must 
be  zero,  thus 


+  V' [VpVl  +  VP  =  0 


(A2 


Expanding  the  first  term  of  Eq.  (A27)  and  using  the  vector 
identity 

V- [VpV]  =  V[7-pV]  +  p(V*V)V 


Eq.  (17)  becomes 

P  +  v{|^  +  v-pvl  +  P(V-V)V  +  VP  =  0  (A28) 

The  second  term  of  Eq.  (A28)  is  zero  by  the  conservation  of 
mass  so  that  Eq.  (A28)  reduces  to 

3\7  —  — 

p  +  p(V*V)V  +  VP  =  0  (A29) 


Using  Eqs .  (A2),  (A5) ,  (All)  and  (A12),  collecting  terms  and 

equating  each  vector  component  to  zero,  Eq.  (A29)  becomes 


(A30) 

(A31) 

(A32) 


Using  the  perfect  gas  relation 


and  the  identity 


IZ  =  p  ^  ‘’-nP 
9  s  9  s 


Eqs .  {A30)  through  (A32)  become 


“  '  ^  *  • 

h-  v'v 

I  4 


M  +  q  = 
9t  ^  9s 


^  9  ^nP 
Y  9s 


(A33) 


U  KJ  ,  _  U 

9t  9s  ~ 


A  9  gnP 
Yq  9n 


(A34) 


ii  +  g  li  =  - 


9t  ^  9s 


A  9  ilnP 

yq  cos  6  9m 


(A3‘9) 


which  are  in  the  desired  form. 


F.  CONSERVATION  OF  ENERGY 

The  conservation  of  energy  for  an  arbitrary  control  volume 


■  •*«  v 


The  rate  of  increase  of 
energy  within  the  con¬ 
trol  volume  with  respect 
to  time 


The  rate  at  which  energy 
is  entering  the  control 
volume  across  the 
boundary 


The  net  rate  of 
work  done  on  the 
fluid  at  the 
boundary 


(A36) 


Neglecting  gravitational  potential  the  energy  per  unit  mass 


where  e  is  the  internal  energy  per  unit  mass  and  q  is  the 
velocity  magnitude.  For  an  inviscid  fluid  in  the  absence  of 
body  forces,  the  only  work  done  on  the  fluid  is  pressure  work 
Thus  Eq.  (A36)  is  written  mathematically  as 


a 

at 


///  [e  +  p  dvol 


-//[e  +^j  pV-dX"  -  //  PV-da  (A3 


Using  Gauss'  Theorem,  Eq.  (A37)  is  written  as 

2  2  _ 
///■^[e  +  ^]  P  dvol  =  -  ///v*  [  (e  +  dvol 

-  Iff  V- [PV]dvol' 


or 


-|^{[e+^]p}  +  V-{[e+^]pv}  =  -V*  [Pvj 

Expanding , 

2  2  2  2 
|^[e  +^]  +  p  |^[e  +^]  +  (PV)V-  [e  +^1  +  [e  +^]  [V-  (pV)  ] 


-V-  [PV] 


L 


2  _  p,  rr2  _  2 

[e  +^]  [||- +V- (PV)  ]  +  p  |^[e  +^]  +  (pV)V*[e+^] 


=  -V*  [PVl 


(A38) 


The  first  term  of  Eq.  (A38)  is  zero  according  to  the  conserva¬ 
tion  of  mass.  Expanding  the  remaining  terms  of  Eq.  (A38) , 


or 


-  i  V.  [PV] 


(A39) 


Equation  (A39)  can  be  written  as 


—  +  —[2—]  = 
Dt  Dt^  2 


-  i(vVP) 


|(V- V) 


(A40) 


With  V  =  qig  in  the  natural  coordinate  system,  Eq.  (A29) 


becomes 


3V  ^  3V 
9t  3s 


-  -  TP 
P 


or 


DV 

Dt 


-  -  VP 
P 


(A41) 
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Taking  the  dot  product  of  v  with  Eq.  (A41) , 


^'Dt  = 


-  ^[V.VP] 


and  substituting  Eq.  (A42)  into  Eq.  (A40) 


Dt  Dt '  2  ' 


v-§|  -  i„.v. 


D(qi^) 


_  /V  Di 

52  i  +  q  —i 
Dt  s  ^  Dt 


_  ^  9i  9i 

Dq  .  ,  ,  s  ,  s  1 

Dt  's 


and  using  Eq.  (A.ll)  and  Eq.  (A12)  , 


Dq?  f96? 

TO  S  ‘J'lt  "n  *  S  at  "n,> 


.  ^2, II  i  t  cos  9  If  i  1 

^  9s  n  9s  m 


Therefore 


-  DV 
V‘Dt 


?  /Dq  ?  ^  .  90  ^  2  901'?  „  rM  ^  i 


5_r3_ 

Dt  '■  2 


Therefore  Eq.  (A43)  reduces  to 


isn 


'  p  — 

:  =  -  ^(V.V) 


From  the  conservation  of  mass 


1^  +  V- (pV)  =  0 


or 


i£ 

at 


+ 


V-  Vp  + 


p(V*  V) 


0 


which,  in  the  natural  coordinate  system  gives 

If  =  -P(VV) 

=  ( 7  •  V) 

Substituting  Eq.  (10)  into  Eq.  (9) , 


or 


at 


i  ££ 

p  Dt 


De  P  Dp 

Dt  ”  2  Dt 

P 


The  definition  of  the  modified  entropy  is 


dS  = 


1 

Y 


do 

T 
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.V.V. 


and  the  first  law  of  thermodynamics  for  an  elemental  mass 
requires  that 


de  +  Pd(i) 


so  that 


dS  = 


-  ;ilde  +  Pd  (i)  1 


=  -  :^!de  +  4  dp) 


and  Eq.  (A46)  implies  that 


=  0 


(A47) 


which  states  that  modified  entropy  is  conserved  along  streamlines. 


G.  MODIFIED  ENTROPY  EVALUATION 


With  the  modified  entropy  defined  as 


-=A 


-  1  '"“r 

V  /  T 


(A48) 


From  the  first  law  of  thermodynamics 


» 


m 

i2i*» 
ifW 
iiS’ 


dQ_  =  de  +  P  dv 
K 


(A49) 


where  e  is  the  internal  energy  and  v  is  the  specific  volume, 
Substituting  Eq.  (A49)  into  Eq.  (A48) , 


=  T  / 


de  +  P  dv 
T 


1  “  ,  i 

-  hi  I  I 


P  dv, 


For  a  perfect  gas 


P  v  =  R„T  and  de  =  C  dT 

Cj  V 


for  which  Eq.  (A50)  becomes 


^^A  - 


T  B  C„  dT  B  R_  dv 

/  -V- "  /  -V-J 

A 


which  after  integration  yields 


Sa  =  Sq  +  iln  +  R  an  — ] 

A  B  Y  '  V  v^-" 


Substituting 


Rg  =  C^(v-l) 


into  Eq .  (A51) 


(A50) 


{A51) 


.•  V  'j-  r.»  r  *  'j-  •>  '>  *>  k>  >>  v- 


V 


H.  FURTHER  TRANSFORMATIONS 


Equations  (A24) ,  (A33)  through  (A35)  and  (A47)  are  in  the 

required  form.  However,  some  of  the  variables  need  to  be 
transformed,  namely,  derivatives  of  ?n  P  need  to  be  expressed 
in  terms  of  A  and  S  and  derivatives  of  p  must  be  expressed  in 
terms  of  q  and  A. 

Using  the  first  law  of  thermodynamics  expressed  as 

dQ  =  dh  -  -  dP 
P 


Eq.  (A48)  may  be  written  as 


-Y  T  VS  =  Vh  -  i  VP 


or  for  the  i^  component 


15.  =  i  ^ 

3s  T  3s  pT  3s 


(A5 


As.suming  to  be  constant  and  using  the  equation  of 
for  a  perfect  gas  Eq.  (A54)  becomes 


state 


15.  =  Ie  l_r_l_i  _  R  ^  P 

3s  T  Ss^R^pJ  ^  3s 


and  since  =  YP/(Y“i) 


3S  _  3  .  P  .  3  iln  P 

-s  (y-1)  P  3s^R„P^  G  3s 
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Since  R  is  a  constant  it  may  be  moved  in  or  out  of  the 
derivatives  giving 


3  r  S  , 
^  9sfR  ^ 

—i—  P  ^  f^j  _ 
Y-1  P  9s^p^ 

9  £n  P 

9s 

and  since  y  is 

also 

constant,  similarly 

3  r  S  1 

= 

y  ^  3  r£li 
v-1  Py  3s  p 

9  in 

3  s 

2 

Using  A  =  yP/P 

3  f  S  , 

-  ■>  9s‘r^I 

Y  1  9A^ 

9  £n  P 

Y-1  ^2  9s 

9  s 

2Ay  3A 

3  £n  P 

(y-.i)A^ 

3s 

and,  on  rearranging 


5  y.n  P 

3  s 


V-I  A  ds  •  9s^R' 


Substituting  Eq-  (ASS)  into  Eq.  (A33)  yields 


M  +  A-[S]  =  0 

9t  ^  ^  9s  ^  Y-1  8s  9s^R^  ^ 


and  the  derivative  of  In  P  has  been  removed. 


since 


p 

logarithmic  differentiation  yields 


For  a  pure  substance 

P  =  P(P,S) 


and  an  isentropic  process 


P  =  P(p)s 


Differentiating 


li)  =  -  .2 

9p^S  dp 


p 


therefore 


P 


Substituting  Eq.  {A58)  into  Eq.  (A57)  gives 


0 


(A67) 


>«  •I.*  !.» W 


as  ^  as 
at  ^  '5  ^  = 


Note  that  Eqs .  (A65)  through  (A67)  are  in  the  form 


ax  ,  ,  ax 
at  as 


z 


If  X  =  X(s, t) ,  then 

ax  ax  , 

dX  =  dt  +  ds 

at  as 


and 


=  i2^  +  1e 

dt  at  dt  as 


(A68) 


where  (ds/dt)  is  the  "characteristic"  direction  in  the  (s,t) 
plane  along  which  Eq.  (A68)  holds.  Along  this  direction ,  the 
equations  may  be  transform.ed  from  partial  to  ordinary  differ¬ 
ential  equations.  What  is  desirable  is  to  find  transformed 
variables  in  terms  of  which  Eqs.  {A63)  and  (A64)  will  be  in 
the  same  characteristic  form.  For  this  reason  the  modified 
Riemann  variables 


Q  =  q  +  As/Rq 


(A69) 


R  =  q  -  AS/R^, 


(A70) 
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are  defined,  where  S/R^^  is  defined  in  Eq.  (A53)  .  The  modified 
Riemann  variables  can  be  introduced  by  manipulation  of  the 
equations  as  follows.  First,  multiplying  Eq.  (A63)  by  S, 
gives 


S  +  -jsl?  +  If  +  ■5''®  «  Ih'  -  " 


(A71 


and  multiplying  Eq.  (A67)  by  A  gives 


A  If  *11  =  0 


(A72 


On  adding  Eqs ,  (A71)  and  (A72)  to  Eq.  (A64),  and  introducing 


A  |a  .  A  |a  .  0 

^  0  S 


(A73 


and 


A  S  1^  -  A  S  14  =  0 

9s  9s 


(A74 


after  appropriate  rearrangement 


9q  ,  9S  _  9A  9q  9S  „  9A  9q 

^+A^  +  S^  +  q^+qA^+qS^+A3f 


at 


at  "  ^  Tt 


3s 


3s 


3s 


-  ^  *511  +  As|A  *  A 


-  -4^  ^  q  A  S  [||-  +  cos  0  = 


la 

3s 


0 


Rearranging  and  introducing  Eq.  (A69)  gives 


It  ^  i  = 


Y-l 


A  [S 


2  ,  .3 


^qA  S  (||+cos  e  ||] 


(A75) 


By  subtracting  Eqs.  (A71)  and  {A72)  from  Eq.  (A64) ,  again  using 
Eq.  {A73)  and  Eq.  (A74) ,  and  introducing  Eq.  (A70) ,  one  also 
obtains 


3R  3R 

at  “I-*’  al 


^>1 


+  =  »  H-’ 


(A76) 


Together  with  Eqs.  (A65)  through  (A67) ,  Eqs.  (A75)  and  (A76) 
form  the  system  of  equations  which  describe  the  isentropic  flow 
of  an  inviscid  perfect  gas  under  unsteady,  compressible  condi¬ 
tions,  and  which  may  be  solved  as  a  system  of  ordinary  differ¬ 
ential  equations  along  characteristic  trajectories  in  the 
space- time  domain. 

I.  SHOCK  JUMP  EQUATION 

The  analytical  expression  for  the  change  in  the  extended 
Riemann  variable,  Q,  across  a  stationary  shock  is  derived 
below . 

Consider  a  shock  moving  with  a  velocity  into  a  fluid 
with  velocity  qp  as  depicted  in  Fig.  A12,  with  the  conditions 


Figure  Al2 


Figure  A13 


to  the  left  of  the  shock  being  denoted  by  A  subscripts  and 
conditions  to  the  right  by  B  subscripts.  If  a  velocity  of  -V^ 
is  imposed  on  the  entire  system,  the  situation  becomes  one  of 
a  stationary  shock,  as  depicted  in  Fig.  A13.  In  both  cases, 
the  high  pressure  side  is  to  the  left  (A  side) .  Since  all 
velocities  are  defined  positive  to  the  right,  a  positive  value 
for  the  relative  mach  number  is  defined  as 


w  =  - 


B 


{A77) 


The  extended  Riemann  variable,  Q,  is  defined  as 


q  +  As 


(A78; 


where 


S 


2 

Y-1 


1 

Y (Y-1) 


p" 


(A79) 


If  all  velocities  are  non-dimensionalized  by  Ap  and  pressures 
and  densities  are  non-dimensionalized  by  their  values  at  B, 


c. 

f.V>‘ 


t 


r: 


'V*  V 


Substituting  Eqs.  (A86)  and  (A81)  through  (A83)  into  Eq.  (A80) 
gives 


APPENDIX  B 


USING  EULER- 1  ON  THE  NPS  VM/CMS  SYSTEM 


A.  MEMORY  REQUIREMENTS 

Storage  should  be  defined  as  one  mega-byte  to  ensure  ade¬ 
quate  memory  for  DISSPLA  requirements.  Also  ensure  that  room 
exists  on  the  user's  disk  for  the  listing  file  if  that 
output  option  is  selected. 


B.  TERMINAL  REQUIREMENTS 

The  type  of  terminal  required  to  run  EULER-1  depends  on 
the  output  option  selected.  When  using  the  tabular  listing 
output,  any  terminal  tied  into  the  VM/CMS  system  may  be 
used.  Graphical  output  may  also  be  created  using  any  terminal 
and  the  graph  stored  on  the  user's  disk  in  a  metafile  for  later 
viewing  at  a  graphics  terminal  or  for  printing.  To  have  the 
graph  stored  on  the  disk,  use  the  "COMPRS"  command  on  line 
223  of  the  program  and  comment  out  the  "TEK618"  command  on 
line  224.  When  graphical  output  is  selected  and  it  is  desired 
to  view  the  graph  at  the  terminal  at  execution  time,  an  IBM 
3277-Tek618  dual  screen  terminal  must  be  used.  To  use  this 
option,  use  the  "TEK618"  command  on  line  224  and  comment  out 
the  "COMPRS"  command  on  line  223.  When  using  this  option,  a 
low  quality  hard  copy  of  the  TEK618  display  may  be  obtained 
using  the  TEK4631  hard  copy  unit  attached. 


‘!ti 


C.  PROCEDURE  AND  COMMANDS 

After  loging  on  to  VM/CMS,  use  the  command  "DEFINE 
STORAGE  IM"  to  increase  the  virtual  memory  as  recommended 
above.  This  will  remove  you  from  the  CMS  environment.  Re¬ 
enter  CMS  with  the  command  "I  CMS." 

Edit  the  user  input  area  of  EULER-1  FORTRAN  as  desired 
for  the  particular  problem  being  run.  When  graphical  output 
is  selected,  comment  out  either  the  "COMPRS"  or  the  "TEK618" 
command  as  desired. 

Compute  the  program  using  the  command  "FORTVS  EULER- 1." 
After  compilation,  an  EULER-1  TEXT  and  an  EULER-1  LISTING 
file  will  reside  on  the  user's  disk.  At  execution,  if  a 
tabular  outpvit  has  been  selected,  it  will  be  sent  to  the 
user's  disk  as  "FILE  FT09F001."  To  give  this  file  a  particu 
lar  name,  say  "EULER-1  LISTING,"  use  the  file  definition 
command 

"F  I  9  DISK  EULER-1  LISTING  A  {  PERM" 

at  compile  time.  An  alternate  and  convenient  means  of  com¬ 
piling  the  program  is  to  set  up  an  exec  file  on  the  user's 
disk  by  the  name  "EULER-1  EXEC"  which  contains  the  commands 

F  I  9  DISK  EULER-1  LISTING  A  (  PERM 

FORTVS  EULER- 1 

Then  to  compile  the  program  and  define  the  output  file,  use 
the  command  " EULER- 1." 

After  compilation,  to  execute  the  program,  use  the 
command  "DTSSPLA  EULER-1."  The  message 


...YOUR  FORTRAN  PROGRAM  IS  NOW  BEING  LOADED... 
...EXECUTION  WILL  SOON  FOLLOW... 

should  appear,  followed  shortly  by 

. . . EXECUTION  BEGINS . . . 

If  tabular  output  was  selected,  proper  termination  will 
result  in  a  ready  message.  If  graphical  output  was  selected 
using  the  TEK618  option,  the  plot  will  appear  on  the  TEK618 
terminal,  the  program  will  pause,  and  a 

"...press  ENTER  to  continue" 

message  will  be  displayed  on  the  3277  terminal.  At  this 
point  the  hard  copy  must  be  made,  if  desired,  using  the 
TEK4631  hard  copy  unit.  After  pressing  the  ENTER  key  on  the 
3277  terminal,  the  plot  will  be  erased  and  the  program  will 
terminate.  Proper  termination  will  be  indicated  by  the  message 


END  OF  DISSPLAY  9.2  ####  VECTORS  IN  1  PLOT... 
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APPENDIX  C 


EULER-1  FORTRAN  CODE 
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